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We report on a numerical study of the density matrix functional introduced by Lieb, Solovej and 
Yngvason for the investigation of heavy atoms in high magnetic fields. This functional describes 
exactly the quantum mechanical ground state of atoms and ions in the limit when the nuclear charge 
Z and the electron number N tend to infinity with N/Z fixed, and the magnetic field B tends to 
infinity in such a way that B/Z A ^ Z — > oo. We have calculated electronic density profiles and ground 
state energies for values of the parameters that prevail on neutron star surfaces and compared them 
with results obtained by other methods. For iron at B — 10 12 G the ground state energy differs by 
less than 2 % from the Hartree-Fock value. We have also studied the maximal negative ionization 
of heavy atoms in this model at various field strengths. In contrast to Thomas-Fermi type theories 
atoms can bind excess negative charge in the density matrix model. For iron at B = 10 12 G the 
maximal excess charge in this model corresponds to about one electron. 
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I. INTRODUCTION 

The properties of matter in magnetic fields of the ex- 
treme strength of 10 12 Gauss and higher have been the 
subject of numerous investigations since the early sev- 
enties, a major impetus being the discovery of pulsars 
in 1968 and the resulting interest in magnetized neutron 
stars. We refer to ]l|, (|], ||, [Q, |]J for general reviews on 
this subject and lists of references. The standard Hamil- 
tonian of atomic physics, 
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is usually taken as a starting point for the study of 
atoms in the atmosphere and outermost crust of neu- 
tron stars. Here N is the number of electrons that 
move in the Coulomb field of a nucleus, localized at 
the origin with charge Ze, and in a homogeneous mag- 
netic field B = (0, 0, B) with vector potential A(r) = 
(1/2) (— yB, xB, 0). The Hamiltonian p) operates on 
antisymmetric wave functions € /\ x L 2 (R 3 ;C 2 ) of 
space and spin variables, and a = (<7i, <j%, (73) is the 
vector of Pauli matrices. Units are chosen such that 
h = e — 2m e — 1, c = 1/a « 137; the energy unit is 



then four times the Rydberg energy, i.e., 54.4 eV. Besides 
the atomic Hamiltonian ([!]) it is, of course, important to 
study the Hamiltonian for molecules and matter in bulk, 
but the present paper is only concerned with ([!]), more 
specifically with its ground state energy 



E Q (N,B,Z)= inf (*,r/ w , BiZ $), 
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where ^0 is a ground state wave function. 

Previous works on matter in strong magnetic fields can 
roughly be divided into two classes. On the one hand 
the focus has been on light atoms, in particular hydro- 
gen with Z = 1, on the other hand on heavy atoms with 
high Z. The present contribution falls into the second 
class. Here Z — 26 plays a special role because iron is 
believed to be the most abundant element in the surface 
layer of a neutron star Q| , [Q . For such heavy atoms it is 
reasonable to expect that important aspects can be ex- 
tracted from an asymptotic analysis in Z, and since 10 12 
G is large even compared with the natural atomic unit 
Bq = m 2 e 3 c/H 3 = 2.35 x 10 9 G, an asymptotic analysis 
in B is equally called for. 1 



1 With our choice of units 2m e = 1 and the magnetic field is 
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The asymptotic behavior of the energy ^ and density 
(^) as iV, Z — ► oo , where N/Z is fixed and the magnetic 
field B is allowed to vary with Z as well, has recently been 
rigorously studied by Lieb, Solovej and Yngvason g], [Q, 
H . In these papers it was proved that the ground state 
properties of (1) can in this limit be evaluated exactly 
by five nonlinear functionals corresponding to different 
physics at different scales of the magnetic field B as mea- 
sured by powers of Z . These five parameter regions are 
characterized as follows: Region 1, B <C Z 4 / 3 ; Region 2, 
B ~ Z 4 / 3 ; Region 3, Z 4 / 3 «B« Z 3 ; Region 4, B ~ Z 3 \ 
Region 5, B > Z 3 . Here B < Z p , B > ZP and B ~ Z'P 
means respectively that the ratio B/Z p tends to 0, oo or 
a constant ^ as as Z — > oo. 

The asymptotic theories corresponding to Regions 1- 
3 are semiclassical theories of Thomas-Fermi type that 
have been extensively applied to neutron stars in the 
past, see, e.g., §, §, §, 0, @, @ and g§. Salient 
features of atoms in region 5 were captured by a different 
density functional theory already in the papers jl4j and 
jL5|. However, the conditions on the surface of a typi- 
cal neutron star correspond rather to region 4, and this 
asymptotic region is also the most interesting one from 
the mathematical point of view. In fact, in Q it is shown 
that it can be described by a functional of a novel type, 
where the variable is not a density, but a function with 
values in density matrices. Moreover, this theory covers 
regions 3 and 5 as limiting cases. We refer to it as the 
density matrix (DM) theory. 

In view of the fact that the DM theory is an exact 
limit of quantum mechanics it is important to know its 
properties in some detail. Being an asymptotic theory it 
is clear that it does not encompass the same information 
as the full Hamiltonian at finite Z and B. In particular 
the DM theory does not capture exchange-correlation ef- 
fects, and it is a theory of very strong fields in the sense 
that all electrons are confined to the lowest Landau band. 
These features should not be considered as a shortcom- 
ing of the DM theory, however. In fact, the hardest part 
of the derivation of the limit theorems in Q is precisely 
to prove rigorously that contributions from exchange and 
higher Landau bands vanish in the limit considered. The 
DM theory should be judged in its own merits: It is enor- 
mously more simple numerically than the full quantum 
mechanical problem (||) and it is a well defined starting 
point for more refined approximations. 

In the present contribution we report on a numerical 
study of the DM theory for atoms. We have computed 
ground state energies and electronic density profiles over 



a wide range of parameters and compared them with re- 
sults obtained by different methods. In particular we 
compare the DM theory to the semiclassical theory that 
applies in Region 3, the simple density functional the- 
ory for Region 5, and also to other density functional 
Q, and HF [jl8| calculations. The difference be- 
tween DM and HF calculations of ground state energies 
is less than 2% where data are available so that com- 
parison can be made. This is remarkable in view of the 
fact that for standard Thomas Fermi theory with B = 
the Scott term, which corrects for the rough treatment of 
the electrons close to the nucleus in TF theory, must be 
incorporated in order to achieve such a good numerical 
agreement, cf. |19 . Thus, at least at this field strength, 
DM theory is closer to HF theory than might have been 
expected. A more precise statement requires an analysis 
of the next to leading order terms in the asymptotic ex- 
pansion of the ground state energy. Such an analysis has 
yet to be carried out. 

Another point where the DM theory differs from semi- 
classical theories is in the possibility of negative ioniza- 
tion. It is a general feature of Thomas- Fermi type the- 
ories, based on potential theoretical arguments, cf. p0[ , 
that the number of bound electrons never exceeds Z. For 
the quantum mechanical problem the meaning of this is 
simply that the binding energy of an excess electron must 
necessarily be of lower order in Z than the ground state 
energy. On the other hand it is known that a magnetic 
field enhances binding, for instance the Hamiltonian (1) 
with N = Z + 1 has infinitely many bound states for 
B / ^T|. In the limit of extremely strong fields in 
Region 5 the negative charge can even be as large as 2Z 
[Q. The only rigorous results on the DM theory con- 
cern this extreme limit, but our numerical computations 
clearly show negative ionization that increases with B. 
It seems, however, that in order to approach the 2Z 
value extremely strong fields are needed; even at fields 
as strong 2 as 10 18 G the excess charge for iron is "only" 
about 23 %. 

Our interest in negative ionization is also motivated by 
its relation to another question, the binding of atoms into 
molecules and chains. Although a rigorous mathemati- 
cal theorem linking these two aspects of binding does 
not seem to exist, it is a fact that in regions 1-3, i.e., 
for B Z 3 , molecular binding energies are vanishingly 
small compared to ground state energies, whereas in re- 
gion 5 binding becomes extremely strong: For a diatomic 
molecule the binding energy is 6 times the ground state 



actually measured in units of 4Bo = 9.40 x 10 9 G. 

2 The computations at these extreme field strength were car- 
ried out mainly to test the mathematical properties of DM 
theory. It is clear that doubts about the applicability of the 
nonrelativistic Hamiltonian (|l|) can be raised in such extreme 
fields, even for very heavy atoms. 
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energy of an individual atom! The question whether iron 
is weakly or strongly bound at field strengths of the order 
10 12 G has been controversial over the past 25 years. The 
best HF results |H indicate weak or no binding, but the 
computations are difficult for they amount to subtract- 
ing one large number from another. It is decisive to treat 
the molecules and the individual atoms consistently by 
the same numerical methods so that unavoidable errors 
cancel as far as possible. Since the DM theory is numeri- 
cally much simpler than HF theory it is easier to achieve 
this in the former and we shall return to the binding 
question in a separate paper. The atomic computations 
presented here are a necessary preparation for the study 
of molecules and chains. 



II. THE DENSITY MATRIX THEORY AND ITS 
LIMITING CASES 

The density matrix theory |(| , Q is based on an energy 
functional, where the variable is a mapping T: rj_ — > T r± 
from rj_ = (x, y) € R 2 into density matrices, i.e., nonneg- 
ative trace class operators on L 2 (R, dz) . In a magnetic 
field of strength B these operators have to satisfy the 
condition 



< T r± < (B/(2n))I 



(4) 



for all rj_. Let T r± (z,z') denote the integral kernel of 
T r± and put Pr(r) = T r± (z,z) for r = (r±,z) G R 3 . 
The density matrix functional for an atom with nuclear 
charge Z is defined by 



'[n = -/ 



d 2 T v± {z,z') 



dz'- 



d 3 r - Z 



Pr{r) 



Pr(r)Pr(r') # # 



d 3 r 



(5) 



In the density matrix theory the electrostatic interac- 
tions are treated classically but the kinetic energy for the 
motion along the magnetic field quantum mechanically 
by the —d 2 /dz 2 term. In directions perpendicular to the 
field the motion is restricted by the "hard core" condition 
(|J). This condition reflects the fact that the density of 
states per unit area for free electrons in the lowest Lan- 
dau band is B/(2ir). The functional (||) is plausible if 
one thinks of T r± (z, z') as an approximation to 

N 

N j ^ (r±,z;r 2 , ...,r N )%(r_ L ,z';r 2 ,...,r N ) JJ dr,-, 

(6) 

where V^o is a normalized ground state wave function. In 
the parameter region B Z 4 / 3 the electrons are con- 
fined to the lowest Landau band, and the Pauli Hamil- 
tonian [(p + A(r)) • a} 2 , restricted to the lowest Landau 
band, is precisely —d 2 /dz 2 . 



The ground state energy for TV electrons in DM theory 

is 

E DM {N, B, Z) = inf{£ DM [r] : / p r {r) d 3 r < N}. (7) 

As shown in Theorem 4.3, there is a unique minimizer 
for this variational problem, i.e., 



E DM {N,B,Z) =£ DM [T D N M BZ ] 



(8) 



with a unique T™ B z . The corresponding density, 
PTb.z, satisfies Jpj^ B ,z = N, if N < N c , and 
JPn^bz = N e, if N > N c , where N c > Z is a number 
depending on Z and B. As explained in the next section, 
the minimization problem (Q) amounts to seeking at each 
r± the lowest eigenvalues and eigenfunctions for a one- 
dimensional Schrodinger Hamiltonian —d 2 /dz + VjP"(z) 
where V™ is the self-consistent potential generated by 
the nucleus and P™ B z - 

The density matrix theory is in fact a two parameter 
theory with parameters A — N/Z and rj = B/Z 3 due to 
the scaling relations 



and 



E DM (N, B, Z) = Z 3 E DM (X, r/, 1) 



PN™ B Ar) = Z i Pl^ 1 {Zr). 



(9) 



(10) 



In particular, the ratio to Z of the maximal number 
of electrons that a nucleus can bind in DM theory, 
A c = N c /Z, is a function of 77 alone. 

The DM theory holds a special position in the study 
of the properties of matter in strong magnetic field be- 
cause it provides an asymptotically exact description of 
the quantum mechanical ground state energy E Q and 
electron density P Q as N, Z and B tend to infinity with 
N/Z fixed and B/Z 4 / 3 — > 00. The following theorems 
are proved in 0], Theorems 1.1 and 8.1: 



Theorem II. 1 Let N, Z 

B/Z 4 / 3 -> 00, then 



co with N/Z fixed. If 



E Q (N, B, Z)/E nM (N, B, Z) 



(11) 



Theorem II.2 Let N, Z and B -> 00 with N/Z = A 
and B/Z 3 = 77 fixed. Then 



Z- 4 P Q N , B ,z(Z-^)^p^ 1 (x) 
in the sense of convergence of distributions. 



(12) 



The shape of atoms in DM theory is discussed in Sec- 
tion IV in connection with Figs. 1 and 2. It should be 
kept in mind that by the limit theorems II. 1 and II. 2 the 
DM theory is a theory of heavy atoms. We have chosen 
iron with Z = 26 as our reference because of its astro- 
physical importance. By the scaling relations (|^) and 
( J10| ) it is simple to transform the results to other values 
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of Z. As seen from the figures the atom is approximately 
spherical when the magnetic held is not too strong (< ca. 
10 11 Gauss for iron), but becomes increasingly elongated 
as the field goes up. In fact, as shown in [Q the limit- 
ing cases 77 — s- and 77 — > 00 of the DM theory can be 
described by simpler theories that we now review briefly, 
referring to and || for details. 

The weak field limit, 77 — > 0, is the Thomas-Fermi the- 
ory for atoms in strong magnetic fields, where only the 
lowest Landau band is taken into account (as in DM the- 
ory). This theory was introduced by Kadomtsev jm and 
studied further in a number of publications, see for a 
list of references. In ||, || it is called the STF theory. 
The density functional is 



£ STF [p] = 



4tt 4 
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p(r) 3 d 3 r- Z 
p(r)p(r') 



d A r 



d 3 rd 3 r' 



\r — r 



(13) 



The precise connection between DM and STF theory is 
given in EL_ Eq. (8.11); if E STF (N,B,Z) denotes the in- 
fimum of (13|) 



with subsidiary condition J p < N , then 



lim E DM (X,r],l)/r] 2/5 



E STF (X, 1, 1). 



(14) 



In STF theory, atoms are spherical with a finite radius 

~ Z- 1 /3( B/Z 4/3 ) -2/5^ 

In the opposite parameter regime, more precisely for 
77 larger than a certain critical value, r/ c , DM theory also 
reduces to a density functional theory. The value of rj c 
depends on A; for A = 1 we find rj c — 0.148, which for 
Z = 26 corresponds to B = 2.44 x 10 13 G. The energy 
functional appropriate for such super strong (SS) fields is 



E ss [p]= J [d^p/dz] 2 d 3 r~Z J 



1 f f P(r)p(r<) 

2 ./ ./ |r - r'l 



d 6 r 



d 3 rd 3 r' 



(15) 



with the subsidiary conditions 



p(r) d 3 r < N 



and 



J p(r) dz < S/(2?r) for all r ± . (16) 
In fact, for 77 > rj c , the minimizer of ([5]) has the form 



(z, z') = ^p™(r ± ,z)^/p™(r ± ,z>), (17) 



and (|) evaluated for T DM is the same as ( |15| ) evaluated 
for p um . Atoms in SS theory have the form of a thin 



cylinder with axis in direction of the magnetic field and 
with a cone-shaped region essentially cut out of its inte- 
rior. The radius is finite, R = \j2ZjB. The extension 
along the field is infinite, but the bulk of the electrons is 
confined within a distance ~ Z^ 1 [ln(B / Z 3 )]" 1 from the 
nucleus. 

An even greater simplification occurs in the extreme 
limit 77 — > 00, which we refer to as the hyperstrong (HS) 
case. In this limit the atom becomes effectively one di- 
mensional and is described by a functional that can be 
minimized in closed form. This functional is 



£ HS [p} = J [d^- p /dz] 2 dz-p(o) 



+ / p(zYdz, (18) 



where p{z) is a one dimensional density and the sub- 
sidiary condition is 



p(z) dz<\ = N/Z. 



(19) 



The connection between the SS and HS theories is as fol- 
lows. Let E SS (N, B 1 Z) denote the minimum of (H|) with 
the subsidiary conditions (|l6|), and let E HS (X) denote the 
minimum of ( |l8| ) with the subsidiary condition (|H)|) . Let 
L(rj) be the solution to the equation 

(77/2) 1 / 2 = L( ? 7)sinh(L( ? 7)/2). (20) 

Then we have 

E SS (N, B, Z) = Z 3 L{r 1 ) 2 E HS (\) + Z 3 0{L(r 1 )). (21) 

There is also a corresponding connection between the 
minimizing densities, p|| B z (r) and Pf s (z) for the two 
theories. Namely, 

[Z'mr 1 I pf, B Ar^[ZL( V )}- 1 z)d 2 r ± ^pf (z) 



(22) 

(in the sense of distributions). 

The function L(rf) behaves like In 77 for large 77, so the 
convergence of E 3S to E ns is rather slow. The main in- 
terest in the HS theory is that Pf s (z) and i? HS (A) can be 
explicitly computed: Writing p^ s as [ipf s ] 2 one has 

= V2(2-A) 
1 ' 4sinh[i(2-A)|z|+c] 

^r^)=V2(2+M)- 1 

where tanhc = (2 — A)/2. Moreover 3 , 

£ HS (A) = -IA+IA 2 ^A 3 



for A < 2 (23) 
for A > 2 (24) 



(25) 



We recall that our energy unit is 54.4 eV 
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for A < 2 and £ HS (A) = £ HS (2) = -1/6 for A > 2. Thus 
A c = 2 in HS theory. Eq. is essentially the statement 
of Theorem 3.5 in j|, but with one refinement: By re- 
placing Inrj in that theorem by L(rj) one obtains a neater 
estimate for the error term. The proof of ( pTf ) and (p2|) , 
which follows closely the proof of Theorem 3.5 in (Jis 
given in the Appendix. 



III. NUMERICAL MINIMIZATION OF THE 
DENSITY MATRIX FUNCTIONAL 

In this section we describe in some detail the numer- 
ical methods used to study the DM theory. The task 
is to minimize numerically the density matrix functional 
(||) under the constraints (Q) and N = J Pr(r) d 3 r. The 
density matrix T r± is trace class and can be expressed 
in the form 

oo 

r ri (z,/)^A^[M^)'< 1 W (26) 

3=1 

where cf>J- L (z) is an orthonormal basis in L 2 (R, dz) for 

each r± and < Aj^ < B/(2n), by condition (§]). For 
the minimizer T DM it turns out that 



A 



B/2n \fj<j^ 
ifj>j c r - 



(27) 



with j < j^ 1 - < oo. In fact, j < j^ 1 - < j c with a j : c < oo 
that is independent of r± (but depends on B), cf. jjj, 
p. 553. The functions satisfy the one dimensional 
Schrodinger equation 



dz 2 



Z 

H 



Pr(r') 



and for each (N, Z, B) there exists a unique jjP 



(28) 

such 
6|) is 



that j < jl± if and only if < /i DM , and 
the solutions of the minimization problem. Our strat- 
egy is to minimize the density matrix functional (^|) by 
iteratively solving the set of nonlinear eigenvalue equa- 
tions (|2^) and determining /i DM such that the constraint 
J Pr = N is satisfied. The eigenvalue equations are in- 
variant with respect to rotation around the z-axis. Hence 
they depend only on |rj_|, but even with this reduction 
( p8| ) yields an infinite number of eigenvalue equations, 
one for each value of |rj_|. We reduce them to a finite 
number by making the |rjj-axis discrete. The DM-atom 
has a finite radius R < Ro — \/2N/B. We therefore 
only have to consider the eigenvalue equations for which 



\r±\ < R- Let N r± be the number of eigenvalue equa- 
tions we choose to work with. Let 



Ro 



N r 



1 



(29) 



We solve ( |28[ ) at the N r± points nAj_, n— 1, .., N r± , on 
the |rj_|-axis. Let 



9± f 1 ifr€((n-l)A_L,nAL] 
n ^ ' ] otherwise 



(30) 



We minimize the density matrix functional (|5|) with den- 
sity matrices of the form 



oo ' j_ 



r r Az',z) = J2 E A f x (*')*^ Ax (*)0n(KI 

j=l n=l 



Let 



dz 2 (n& ± ) 2 + z 2 ' 



(31) 



(32) 



and let -0™ denote the eigenfunctions of h n and /t™ the 
corresponding eigenvalues, so that h n ip™ = fi^ip™. We ex- 
press (j> 1 j A± i z ) m terms of approximate eigenfunctions of 
h n which correspond to the Nf, lowest eigenvalues. Hence 
we write 



(33) 



where ip™{z) is an approximation for We deter- 

mine approximations for the basis functions ipf and their 
eigenvalues fi™ by the method of finite elements (FEM) 
for eigenvalue problems, dividing the interval [— z m ,z m ], 
z m > into elements 4 and choosing a polynomial ba- 
sis of degree 5 within each, cf. J22|. Let these solutions 
be ipf(z) where we induce the boundary condition that 
V>™(±z m ) = 0. The eigenvalue corresponding to ^"(z) 
is denoted by . Let N z be the number of samples of 
i/;"(z) values we choose to work with along the z-axis and 
define 



A z = 



1 



(34) 



Then the samples we work with are 4>™{{l — 1/2)A Z — z m ), 
I = 1, . . . , N z . We use the values /x™ as approximations 
for /t™. This basis is chosen because it is expected to be 



close to the solutions 



(z) and as such is a natural 



4 The number of elements we use are 50 — 60, and their width 
varies such that the smallest elements are closest to the origin. 
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starting point for the self-consistent iterations. We have 
now defined the set over which we numerically minimize 



We solve the set of nonlinear equations (E2q) in a self- 
consistent manner iteratively. We define the chain of 
potentials 



V± k) {r) = (1 -a)V^ 1] {r) + aV ( * -1) (r) (35) 



r(fc-l)/ 



for k > and ^ (0) = 0, a £ (0, l] 5 . Let 
the eigenfunctions of the operator 

H nA ± ,(k) = hn + y(k)^ 

with 



jnAi,(k) 



V<; k \r) = 



dV 



where 



L<«(z',z) = 

^«Ai,(fc)^nA x ,(fc)^*^nA_L,(fc) 



(z) be 
(36) 

(37) 
(38) 



Note that 



i»Ai,(0), 



(z) = tpj'(z). With an appropri- 
ate choice of a (we use a € [0.01,0.1]) the sequence 
i^" Ai '' (?) turns out to be convergent and 



lim 

/c — >oo 



b j >{z) = cj) j ± (z) 



(39) 



dV = 



(42) 



/•CO 

Jo 



r 2 + r' 2 
2rr' 



where Q„_ i is an associated Legendre function. The di- 
rect integration of (p7|) is very slow but it is on the other 

(k) 

hand very accurate. To determine Vq at interior grid 
points faster numerically we do the following: We note 
that (37) is the solution of the Poisson equation 



(k) 



r 



(fe) 



(43) 



With the boundary values determined by the direct in- 
tegration we use standard five point difference approxi- 
mation to the Laplacian in order to determine at 
the interior points. We solve the finite difference scheme 
by simultaneous over-relaxation. We find that an over- 
relaxation coefficient of 1.8 yields fast convergence for 
this problem for our choice of Nr ± and N z . When we 

use the finite difference scheme we predetermine Vq 
at a few interior grid points by direct integration. We 
then perform the over-relaxation iterations until we ob- 
tain agreement with the predetermined values up to a 
desired accuracy. We choose to iterate until the first 5 
digits of the solution are identical to all the predeter- 
mined values. 

Now we are ready to go through one iteration of the 
self-consistent iteration scheme in detail. Let us consider 
the A:-th step in the iteration. At the start of this step 



we know F 



(fc-i) 



(recall that 



,nA x ,(0) 



ipf). First we 



J T 3 ■ 

, , and then Vf according to the scheme 



determine V^ k 1 ' 



in L 2 ([—z rn ,z m ]). This defines our self-consistent itera- described above. We next determine the matrix elements 



(|38|) we work with the potentials and of the 



tion scheme. To be consistent with the discrete form of 
form 



I in 



>(*) _ 



5i m ^+Urv} k) (r)^dz (44) 



V(z,r x ) ^^Wj'^nW 



1 = 1 71=1 



where 



1 if 2 6 ((Z-l)A z -z m ,ZA 2 
otherwise 



(40) 



(41) 



for each n = 1 , . . . , N± . Now the eigenvectors of the 
matrices i?^f correspond to the coefficients & in 



(45) 



To determine ( |37| ) we calculate the boundary values on 
the Nr ± x N z grid we work on by direct integration, 



and we denote their corresponding eigenvalues by 



£ ; A - (fe) . I.e., 



5 Self-consistent iterations of the kind discussed here are in 
general not convergent. The coefficient (1 — a) acts as a damp- 
ing factor on the iterations. The value of a is in practice 
chosen by trial and error as high as possible without induc- 
ing instability. That way the iterations converge as fast as is 
possible. 
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H nA ^ k U] lA±Ak) {z) = e ; Ai ' (fc V; Ai,(fc) (^)- (46) 

To determine the eigenvectors and eigenvalues of the 
matrices -ff^ we use the eigen-routines from [ p3| . 
Finally, at the end of the iteration step we determine 
^DM.ffc) sudl that Jr^ ) ± (z,z)d 3 r = N. We are then 
ready for the next step of the iteration. 

We continue the iterations until £ DM (T^) "stops" 
changing. More precisely we choose to stop when the 
change in £ DM (r( fc )) between iterations is in the 7-th 
digit. 

We calculate with N ± = 101, N z = 201 - 301 and 
Nb = 30 — 60. We chose z m in such a way that 2z m is at 
least three times longer than the length of the atom the 
calculation yields. We obtained this criterion by increas- 
ing z m until the ground state energy we obtained became 
stable in the first 6 digits. 

To calculate the energy £ DM we note that 

3 

N b N r ± 

3 = 1 n=l 

= K mi - ,4™ + 2i? DM . (47) 



E UM = 1(0 + E") (50) 
r um = E „ + a d„ (51 ^ 

R DM = \{E' - E"). (52) 

This is how we evaluate the ground state energy and the 
terms it is composed of. Regarding the accuracy in the 
evaluation of K DM , one should be aware of the fact that 
A DM is in general a lot larger than K DM . In the case of 
the STF-theory A STF = 15K STF which is the low mag- 
netic field strength limit for the DM-theory. Therefore a 
slight relative error in A DM yields a much larger relative 
error in K DM . Based on the information given above, we 
estimate the numerical error of the scheme to be about 
±1 in the fourth digit of the ground state energy. We 
therefore show the first 4 digits when we present our re- 
sults of the energy. 

The estimate of A c is done in the following way. 
E UM (A, 77) is calculated as a function of A which is a 
strictly convex function in our approximation, due to the 
finite box. The minimum of this function then deter- 
mines A c . However, the function is extremely flat around 
the minimum so it is difficult to determine its position. 
We make a linear approximation of -E DM (A, rf) using two 
close lying points below the minimum. The value thus 
obtained gives a lower bound on A c close to the true value. 



Here K nM , A° M and i? DM are respectively the kinetic, at- 
tractive and repulsive parts of £ DM . With our choice of 
basis -0.™ and V(r) as the attractive potential we have 



E" := 



Since 



d 2 IYjz,z') 
dz' 2 



E / WW 



+ V{r)T r ^z,z) 
d 2 



d A r 



dz 2 ' ^ 



<t>i ± (z) d 6 r 



N b N b N b N r ± 
j=l 1=1 fc=l n=l 



N b N b N b N r ± 
j=l 1=1 fc=l n = l 

N b N b N r ± 

= EEE 2 ^ A i A " Ax (4)V 

j = l 1=1 71 = 1 



.4 r 



V(r)T r± (z,z) d 3 r 



(48) 



(49) 



we obtain 



IV. ATOMIC PROPERTIES IN DM THEORY 

The results of our numerical computations are pre- 
sented in Tables I- VI and illustrated in Figures 1-5. As 
remarked before, iron is of special importance in astro- 
physical context, and for this reason we state our results 
for the reference value Z = 26. The scaling relations (||) 
and (|lo| ) allow an easy transformation to other values. 
Moreover, in this section magnetic fields are measured in 
Gauss and energies in ke V to facilitate comparison with 
astrophysical data and other computations. To trans- 
form into the units in which the original Hamiltonian (Q) 
is written it should be kept in mind that there the energy 
unit is 54.4 eV and the unit for magnetic field strength 
is 9.40 x 10 9 G. The dimensionless parameter n — B/Z 3 
has for Z — 26 and B = 10 12 G the value 6.053 x 10~ 3 . 
This may seem small for Region 4, but as discussed in [BJ , 
the relevant semiclassical parameter is really 77 1 / 5 and at 
the quoted values of Z and B we have 77 1 / 5 = 0.36. 

Figures 1 and 2 show contour plots of the electronic 
densities of iron atoms with N = Z according to DM 
theory at four different field strengths, ranging from 10 11 
to 10 14 Gauss. It is evident that in the weakest field the 
bulk of the electrons is spherically distributed around 
the nucleus. With increasing field strength the spheric- 
ity gets more and more distorted. The atom is composed 
of cylindrical shells that decrease in number as the field 
goes up. Between 10 13 and 10 14 Gauss a transition to a 
cylindrical shape with a single shell takes place. The 
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number of shells corresponds to the number of eigen- 
functions of the one dimensional Schrodinger operators 
—d 2 /dz 2 + V™(z) that contribute to the density ma- 
trix in the sum (|26|). The critical value, r) c , at which 
this number has dropped to one, was determined numer- 
ically for A = 1 to be rj c = 0.148. This corresponds to 
B = 2.44 x 10 13 G for Z = 26. 

In Table I the ground state energy ^ of iron is shown 
as a function of the field strength for B between 10 10 and 
10 14 Gauss and for various values of the ratio A = N/Z of 
electron number to nuclear charge. A comparison of some 
of these values with results obtained by other methods is 
given in the next section. 

Table II shows the results tor A c as a function of B in 
DM theory. At the extremely strong field of 10 18 Gauss, 
corresponding to 77 1 / 5 = 5.7, one finds A c = 1.232, which 
is still quite far from the HS value A c — 2. However, 
compared with Thomas Fermi theories, where A c is al- 
ways 1 p0| , the negative ionization is noticeable even 
for the weaker fields. The value A c = 1.046 for iron at 
B = 10 13 G corresponds to an excess negative charge of 
1.2 electrons. The binding energy of the excess charge in 
DM theory is shown in Table III. 

As remarked at the end of Section III a precise deter- 
mination of A c is difficult and the values quoted should be 
regarded as lower bounds at the respective field strengths. 



V. COMPARISON WITH OTHER THEORIES 

As discussed in Section II, DM theory simplifies in the 
two limits, 77 — s- and rj — > 00, which are respectively 
described by the STF and HS theories. In order to study 
the rate of this convergence we have compared the ground 
state energies at A = 1 in these three models in Table 
IV and plotted them in Fig. 3. It is remarkable how 
closely the STF ground state energy approximates the 
DM energy even at fields as strong as 10 15 Gauss, while 
the electronic densities in DM theory deviate appreciably 
from the spherical shape of STF theory already at 10 12 
Gauss as seen in Fig. 1. Thus in this case at least, the 
energy calculations are much less sensitive to the details 
of the model than density calculations. There is, how- 
ever, another way of comparing the densities in STF and 
DM theories. In Fig. 4 we have plotted together the STF 
density and the spherically averaged DM density, and one 
sees that they are quite close for the bulk of the electrons, 
up to fields of the order 10 12 G. 

It is apparent from Fig 3 that the DM energy values 
approach the HS values as the field goes up, but the con- 
vergence is very slow and the asymptotic regime has not 
yet been reached at the strongest fields considered for the 
DM computations. On the other hand, it is interesting 
how well the HS density (|2^) fits the DM density inte- 
grated over the cross section of the atom as shown in Fig. 
5. 

Finally, in Tables V and VI we compare the ground 



state energy E computed in the DM and STF theo- 
ries with values obtained by some other methods in the 
literature. The comparison is made for iron at 10 12 
G, since this case has been considered in a number of 
sources. In table V the value for DM theory is com- 
pared with Hartree Fock (HF) flql , density functional 
IH , (l7) , restricted variational (RV) |3L Thomas-Fermi- 
Dirac (TFD) |]| and STF calculations]!]. In Table 5 the 
splitting of the ground state energy into its various parts 
(kinetic, attractive, repulsive, exchange) is compared for 
DM, STF and HF calculations. 



VI. CONCLUSIONS 

We have carried out a numerical study of the den- 
sity matrix model that describes exactly the quantum 
mechanical ground state of atoms in a homogeneous 
magnetic field in the asymptotic limit when the nuclear 
charge Z, the electron number N and the magnetic field 
B tend to 00 with N/Z fixed and B/Z 4 ' 3 -> 00. The 
calculations demonstrate the following features of heavy 
atoms in high magnetic fields as the field strength in- 
creases: A transition from an approximately spherical 
shape to a highly elongated shape, accompanied by a de- 
crease in ground state energy and increasing ability to 
bind excess electrons. We have also compared the DM 
model with the semiclassical Thomas-Fermi theory and 
the one dimensional density functional theory that de- 
scribe respectively its low and high field limits. When 
(B/Z 3 ) 1 / 5 is of order unity these simpler theories are nu- 
merically and conceptually wrong and the full DM theory 
should be used. The quantitative agreement of DM the- 
ory with Hartree Fock calculations is quite good in strong 
fields; for iron at B = 10 12 G the difference in binding 
energies is less than 2 %, The DM theory, however, is 
much simpler computationally than HF theory for large 
atoms and appears suitable as a starting point for more 
refined approximations and for the study of molecular 
binding in strong fields. 



APPENDIX: 

We give here a proof of Eqs. ( |2l| ) and (p2|). It con- 
sists essentially in an improvement of the estimates in 
Proposition 3.3 in pj. 

With L(rf) defined by ( pp| ) we define for any density p 
a rescaled density p v by 

p(r x ,z) = Z 4 V L( V )p v (Z V ^ 2 r ±1 ZL( V )z). (Al) 

We can then write the SS functional ( fl5l) as 

£ ss [p] = Z 3 L( v ) 2 £f 1 s [p 11 ], (A2) 

where £^ s is defined by 
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£ T [Pv\ 



dz 



d A r 



p v (r)V v (r)d 3 r (A3) 



p„(r)y,(r-r / )/>,(r / )d 3 rd s r / 
with the rescaled Coulomb potential 

V n (r) = Lfa)- V + r?- 1 Lfa) a ri) 1/2 - (A4) 
The functional £^ s has a ground state energy 

E s v s (X)=inf{^[ Pri ]:p v eC ss ,J p v < A, 

p„(rj.,«)d»<l}. (A5) 



and a corresponding minimizing density denoted by by 
p^ s . It is related to p ss by the scalin g (|Al| ). The energy 
E^ s is related to E ss by the scaling flA2|), 



E SS (N, Z, B) = Z L(r)) E^ s (\) 



(A6) 



The improvement of proposition 3.3 in [Q] is stated in the 
following lemma: 

Lemma A.l For any choice of X and T there is a con- 
stant C(A, T) such that 



J p(r x ,0)d 2 r x - j V vP 



<C{\,T)L{rj)- x (A7) 



holds, provided p > satisfies J p < A, p(r±, z) = for 
\r±\ > and 



We can choose C(A,T) 
8y/nT^ 2 ln(V2 + \/6). 



d 3 r < T. 



A + Sv^A 1 / 4 ^/ 4 



Proof : Following the proof of proposition 3.3 we 
write the difference on the left side of ( A7) as A1+A2+A3 
with 



(A8) 



* = -/ ^(rWrJd'r, 

/|7\l|>1 



A 2 = / K,(r)[p(r ± ,0)-p(r)]d 3 r (A9) 
'|rx|<i 



and 



A, 



1- / V„(r)dz 

|7\l|<1 



p(r±,0)d 2 r ± . (A10) 



With the same arguments as in g we obtain 
1*1 < ' 



L( v y 



and 6 
1*1 



V2 



1*1 < ^a 1 / 4 ^ 4 



S inh- V/ 2 /(£MKD) - 1 



(All) 
(A12) 

(A13) 



p(r_L,0) d 2 r ± 



s\nh~ l {ri x/2 /{L{ri)r)) - 1 



2 V2 



Now we deviate from Q. Estimating the integral in 
(|A13j ) we obtain 



nl/2r,l/4 



|A 3 | < 2V2ttT l/z 2 
x sup 

0<r<%/2 



Mv) 

< 8V¥T 1/2 ln(V2 + V6)/L(n). 



sinir 1 (f]^ 2 /(L(r])r)) - 1 



,1/2 



(A14) 



The last inequality comes from the following. If a = 
V^/ L (v), then 

2sintT 1 (a/\/2) - L{rj) = (A15) 

by the definition @ of £(77). Using ( |A15| ) we get 

(2sinh" 1 (a/r) - L(n))r 1/2 



sup 
sup 

0<r<V2 



2sinh~ 1 (a/r) - 2 sinh _1 (a/\/2) 



2sinh _i (a/V2)-L(77) 



.1/2 



= sup 

0<r<\/2 

= 2 sup 

0<t-<V5 

< 2 sup 

0<r<V2 



2sinh~ 1 (a/r) - 2 sinh _i (a/V2 



,1/2 



ln(V2/r) + In 



ln(V / 2/r)r 1/2 



+ \/a 2 + r 2 



+ Va 2 + 2 



a/2 



+2 sup 

0<r<^2 



In 



a + \Jo? + 2 



< 2 5/4 In 2 + 2 5/4 ln((l + Vs) /V2) 
= 2 5 / 4 ln(v / 2 + \/6). 



(A16) 



In y|, p. 541 there is a power of 1/2 missing on T in the 
estimate for A3. 
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This concludes the proof. 

In analogy with proposition 3.4 in [Q the following 
lemma is a corollary of Lemma A.l: 

Lemma A. 2 For p as in lemma 1 there is a constant 
C (A,T) such that 

p(r)V v (r - r')p(r') d 3 r d 3 r' - J p{zf dz 
^CWT)^)- 1 (A17) 
where p(z) = J p(r±,z)d 2 r±. 

The proof of (^l|) and ( p2f ) is now identical to the proof 
of theorem 3.4 in B with lemmas 1 and 2 in place of 
Propositions 3.3 and 3.4. 

In order to illustrate the difference between the func- 
tion L(ij) and the approximation L(n) In 77 used in 
the ratio £(77)/ In 77 is plotted in Figure 6. 
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B (G) 

FIG. 1. Contour plots of the electronic density of iron FIG. 3. The ground state energy of iron atoms as a func- 

atoms in DM theory for B = 10 11 Gauss (left) and B — 10 12 tion of the magnetic field strength B in DM theory (crosses), 
G (right). The outermost contour encloses 99 % of the neg- STF theory (short dashes) and HS theory (long dashes), 
ative charge, the next 90 %, then 80 % etc., and the two 
innermost 5 % and 1% respectively. 
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FIG. 2. Contour plots of the electron density of iron atoms 
in DM theory for B = 10 13 Gauss (left) and B = 10 14 G 
(right). The contours are drawn in the same way as in fig. 1. 
At B = 10 14 the DM model has simplified and the density is 
described by the SS functional (12). 
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FIG. 4. Comparison of the electron density in STF theory 
(dashed curve) and the spherically averaged density in DM 
theory (solid curve) for B = 10 11 Gauss (right) and B = 10 12 
G (left). 
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FIG. 5. Comparison of the electron density in HS theory 
(dashed curve) and the integral over ri of the density in DM 
theory (solid curve) for B = 10 13 Gauss (a) and B = 10 18 G 
(b). 
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FIG. 6. The ratio L(t])/lnr]. 
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TABLE I. Ground state energy (in keV) of iron atoms 
(Z=26) as a function of the magnetic field B (in G) and the 
ratio A = N/Z of electron number to nuclear charge. 



A\ B 


10 10 


10 11 


10 12 


10 13 


10 14 


0.1 


-3.202 


-8.188 


-20.80 


-52.58 


-122.4 


0.2 


-4.753 


-12.07 


-30.56 


-77.56 


-185.2 


0.3 


-5.838 


-14.80 


-37.47 


-95.04 


-230.2 


0.4 


-6.670 


-16.88 


-42.64 


108.0 


-264.0 


0.5 


-7.277 


-18.41 


-46.56 


117.9 


-289.4 


0.6 


-7.744 


-19.58 


-49.53 


125.2 


-309.3 


0.7 


-8.102 


-20.42 


-51.73 


130.6 


-323.7 


0.8 


-8.321 


-21.01 


-53.12 


134.1 


-333.5 


0.9 


-8.475 


-21.41 


-53.85 


136.6 


-339.9 


1.0 


-8.521 


-21.47 


-54.38 


137.8 


-342.7 


TABLE II. The ratio A c : 


= N c /Z of the maximal negative 


charge 


to nuclear charge as 


a function of B (in 


G) for iron 


(Z=26) 












B 




A c 


B 




A c 


10 1U 




1.020 


10 ib 




1.110 


10 11 




1.026 


10 16 




1.153 


10 12 




1.035 


10 17 


1.184 


10 13 




1.043 


10 18 




1.232 


10 14 




1.061 








TABLE III. The binding 


energy at maximal negative ion- 


ization in DM theory. 








B 


E DM (X C )-E DM (1) 


B 




(A C )-E DM (1) 




E DM (1) 






E™(1) 


10 iu 




0.0011 


10 ib 




0.0073 


10 11 




0.0008 


10 16 




0.0123 


10 12 




0.0035 


10 17 




0.0232 


10 13 




0.0038 


10 18 




0.0203 


10 14 




0.0050 









TABLE IV. Comparison of the ground state energy (in 
keV) of iron atoms in DM theory, STF theory and HF theory 
at various field strengths B (in G). See also Fig. 3. 



j^DM E STF E HS 

W 1 -21.47 -21.57 

10 12 -54.37 -54.07 

10 13 -137.8 -135.8 

10 14 -342.7 -341.2 

10 15 -786.3 -857.0 -0.3346 

10 16 -1623 -2153 -535.3 

10 17 -3065 -5405 -1797 

10 18 -5264 -13583 -3913 



TABLE V. Ground state energy (in keV) of iron atoms at 
B = 10 12 G according to DM theory, HF theory [17], DF com- 
putations, denoted DF a [15] and DF b [16], RV computations 
[23], TFD theory [12], and STF theory [8]. 



DM 


HF 


DF a 


DF 6 


RV 


TFD 


STF 


E -54.38 


-55.10 


-56.10 


-58.3 


-53.13 


-56.21 


-54.07 



TABLE VI. The composition of the ground state energy E 
(in keV) of iron atoms at B = 10 12 G in STF theory, HF the- 
ory and DM theory. K is the kinetic energy, A the attractive 
potential energy due to the nucleus, R the energy of Coulomb 
repulsion and _E ex the exchange energy. 





E 


K 


A 


R 




STF 


-54.07 


10.81 


-97.33 


32.44 





HF 


-55.10 


10.6 


-95.4 


32.7 


-3.06 


DM 


-54.38 


10.43 


-96.90 


32.09 
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